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Mounting evidence suggests that the TeV-PeV neutrino flux detected by the IceCube telescope 
has mainly an extragalactic origin. If such neutrinos are primarily produced by a single class 
of astrophysical sources via hadronuclear {pp) interactions, a similar flux of gamma-ray photons is 
expected. For the first time, we employ tomographic constraints to pinpoint the origin of the IceCube 
neutrino events by analyzing recent measurements of the cross correlation between the distribution 
of GeV gamma rays, detected by the Fermi satellite, and several galaxy catalogs in different redshift 
ranges. We find that the corresponding bounds on the neutrino luminosity density are up to one 
order of magnitude tighter than those obtained by using only the spectrum of the gamma-ray 
background, especially for sources with mild redshift evolution. In particular, our method excludes 
any hadronuclear source with a spectrum softer than E~^'^ as a main component of the neutrino 
background, if its evolution is slower than (1 -|- z)^. Starburst galaxies, if able to accelerate and 
confine cosmic rays efficiently, satisfy both spectral and tomographic constraints. 

PACS numbers: 95.85.Pw, 95.85.Ry, 98.70.Rz, 98.70.Vc 


Introduction. —The discovery of the PeV neutrinos by 
IceCube mm has launched the era of high-energy neu¬ 
trino astronomy. The current data set is compatible with 
a flux in excess with respect to the atmospheric back¬ 
ground, with an isotropic allocation of events on the ce¬ 
lestial sphere and flavor equipartition [Ml- Due to the 
current low statistics, the origin of the high-energy Ice- 
Cube events is not yet known, but an extragalactic and 
mostly diffuse origin appears to be favored BIS]. 

The high-energy neutrino production from cosmic ac¬ 
celerators has been subject of a cascade of theoreti¬ 
cal studies, especially after the IceCube results were 
announced m iH]- Many papers discuss the neutrino 
emission from one specific source class by adopting a 
model-dependent approach, for active galactic nuclei 
(AGNs) mn], star-forming galaxies [TM2H] , gamma-ray 
bursts [2M36] . galaxy clusters [dTHdf)] . and dark matter 
decays [11M5] . 

Alternatively, a more generic approach focuses on the 
phenomenological aspects of the potential sources. For 
example, assuming photomeson production (py) of neu¬ 
trinos, Ref. [15] obtained constraints on the source size 
and magnetic field strength needed to match the Ice- 
Cube flux. Reference [17] hypothesized that the TeV- 
PeV neutrinos were generated via hadronuclear interac¬ 
tions {pp) and concluded that the cosmic ray spectrum 
of the dominant neutrino sources should be harder than 
This is because the associated gamma-ray spec¬ 
trum will extend down to GeV energies, where the flux of 
the isotropic gamma-ray background (IGRB) measured 
with the Fermi Large Area Telescope (LAT) [l8] cannot 
be overshot. The connection with sources of ultrahigh- 
energy cosmic rays has also been considered [iSHST] . 

In this Letter, we complement the existing model- 
independent investigations of pp neutrino sources by 
proposing an entirely new method: Tomographic con¬ 
straints, up to now adopted in studying IGRB sources. 


We base this approach on the measurements of Ref. [52] , 
which analyzed the IGRB data and found that they were 
spatially correlated with galaxy distributions. Gompared 
to the commonly adopted spectral analysis, the tomo¬ 
graphic method allows to efficiently extract a dominant 
IGRB component in certain redshift ranges following 
galaxy catalogs, as originally proposed for dark matter 
detection [581455] . This provides stringent constraints on 
astrophysical sources [52] and dark matter [55] . 

We show that the tomographic approach allows to 
tightly constrain the redshift evolution and the energy 
spectrum of any class of astrophysical source producing 
high-energy neutrinos through pp interactions, especially 
if the source luminosity density mildly evolves as a func¬ 
tion of redshifts. It provides constraints on the expected 
neutrino flux that are more stringent by up to one or¬ 
der of magnitude with respect to the common spectral 
approach (e.g.. Ref. m)- We find that any source with 
a spectrum softer than E~^-^ is excluded, if its redshift 
evolution is slower than (I -|- z)^. On the other hand, 
sources with hard spectrum and fast evolution can still 
be dominant in both gamma-rays and neutrinos. 

Besides the pp origin of the high-energy neutrinos, we 
assume that: (i) the energy spectrum is a power law, 
E~°‘, extending up to PeV energies; (ii) the source lumi¬ 
nosity density evolves as (I -I- up to ^c, and is constant 
for z > Zc, and (iii) the astrophysical sources trace the 
underlying dark matter distribution. The third assump¬ 
tion is generic for any known extragalactic source that is 
likely associated with large cosmic structures. We adopt 
cosmological parameters from Ref. m- 

Gamma-ray intensity. —We first introduce a differ¬ 
ential gamma-ray intensity, Ij{E), as the number of 
gamma-ray photons received per unit area, unit time, 
unit solid angle, and unit energy. It is computed as an 
integral of the gamma-ray window function, Wj{E,z), 
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over the comoving distance x' 

I^iE) = J dxW^{E,z) , (1) 

4nAE'i,AEL.n) (1 + -)“ 
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where n{z) is the source number density at z, {Lj(z)) is 
the mean gamma-ray luminosity emitted between = 
0.1 GeV and i^^ax = 100 GeV in the source rest frame 
(as represented by '), and 


A = 


1 (-E„ax^g„in) - £qj. a ^2 , 

ln(i^max/^min) for ^ = 2 . 


(3) 


The source luminosity density is assumed to evolve as 


n{z){L^{z)) = X 


(1 -I- zY for z < Zc , 
(1 -I- ZcY for z > Zc . 


(4) 


The constant evolution above z^ is motivated by the ob¬ 
servations of infrared luminosity density of star-forming 
galaxies (e.g., [58]). We note that unless the redshift 
dependence continues to increase steeply up to high z, 
our conclusions are largely unaffected. Very-high-energy 
gamma rays are subject to absorption by the extragalac- 
tic background light (EBL). This is taken into account 
through the exponential term in Eq. ([^, where t{E,z) 
is the optical depth [59]. 

For each set of (a, S, Zc), by taking £j q as a free pa¬ 
rameter, we compute the x^ statistic as follows: 

^2 _ I dat 

i ^ 

where /i,dat and (Ji.dat are the spectral intensity data and 
the associated root-mean-square error in the Tth energy 
bin, respectively, and /qth('^^ 7 ,o) is the theoretical model 
intensity for £^fi. The 95% confidence level (CL) upper 
limit on £^ q is obtained by solving — Xmin = 

2.71. 

The top panel of Fig.[2shows the gamma-ray spectrum 
for a = 2.2, (5 = 2, and Zc = 1.5 (blue dotted), compared 
with the IGRB measured by Fermi [48] . The value of the 
local luminosity density ^..y.o corresponding to the 95% 
CL upper limit is = 2-5 x 10^^ erg yr"! Mpc"!. 

Cross correlation with galaxy catalogs .—The cross¬ 
correlation angular power spectrum, G/®, between the 
gamma-ray intensity, and the galaxy surface den¬ 

sity, Eg(n), is related to the angular correlation function 
through the following relation (e.g., [53]): 


,th (^7,0 IO:, Sj Zc) 
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{6E{n) SEJn + e)) =Y PYcos0) , (6) 
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where SI^ = — (I^), = Eg — (Eg), Pe{cos6) is 

the Legendre polynomial, and Wi is the beam window 




FIG. 1. Top: Gamma-ray (blue) and neutrino (red) intensi¬ 
ties for a model with a = 2.2, 5 = 2, and Zc = 1-5. The dotted 
curves correspond to the 95% GL upper limit due to the Fermi 
spectrum data (the green band represents the systematic un¬ 
certainty due to the subtraction of the Galactic emission |48|1. 
The solid curves correspond to the same limit but due to 
the cross-correlation data. IceGube data for the neutrino in¬ 
tensity are shown above 10 TeV, whereas the orange band 
represents the 68% GL region of the corresponding best-fit 
single power-law model [5]. Bottom: Gross-correlation angu¬ 
lar power spectrum between the Fermi data, above 1 GeV, 
and the 2MASS galaxies, compared with the measurements 
by Ref. [52]. Model parameters as well as line types are the 
same as the top panel. 


function (i.e., the Legendre transform of the point spread 
function of the Fermi-LAT [52]). 

The angular cross-power spectrum G^® is computed as 

(e.g., [53]) 

= J ^W,iz)W,{z)P,, (^k =A) , (7) 

where 14% (z) is the integrated gamma-ray window func¬ 
tion, and ITg(z) is the galaxy window function that 
is related to the galaxy redshift distribution, dNg/dz, 
via Wg(z) = {d\nNg/dz){dz/dx)■ We approximate the 
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cross-correlation power spectrum between the gamma- 
ray emitters and the galaxy catalogs as P^g « b^bgP^, 
where Pm is the nonlinear matter power spectrum com¬ 
puted with the publicly available CLASS code [SD], and 
bg and b^ are the bias factors for the catalog galaxies 
and the gamma-ray emitters, respectively. We assume 
that gamma-ray sources are unbiased tracers of the dark 
matter distribution, i.e., b^ = 1. Since astrophysical 
sources are typically positively biased dark matter trac¬ 
ers (e.g., and references therein), it is a conservative 
assumption. 

The cross-correlation analysis of Ref. [55] adopted five 
different catalogs: Two Micron All Sky Survey (2MASS), 
quasars in the Sloan Digital Sky Survey (SDSS), the 
SDSS main galaxy sample, luminous red galaxies in 
SDSS, and radio galaxies in the NRAO VLA Sky Sur¬ 
vey (NVSS). Each of these catalogs traces underlying 
dark matter distribution in a certain redshift range with 
a characteristic bias bg as in Ref. [52]. Although some of 
them represent AGNs, they can be used the same way as 
galaxies, for which we call them “galaxy” catalogs col¬ 
lectively. We use a redshift distribution dNg/dz and a 
typical bias bg appropriate for each catalog, and three 
different energy ranges for the gamma rays (> 500 MeV, 
> 1 GeV, and > 10 GeV) [52|. 

Similarly to the spectral analysis, for each given set of 
{a, S, Zc), we compute the follows: 

*“ = E E - Th‘), (CI‘ - . 

7,g 

( 8 ) 

where 7 and g run through three energy bins and five 
galaxy catalogs, respectively, I and I' represent the mul¬ 
tipole bins of the measurements, and Gov is the covari¬ 
ance matrix. We again use = 2.71 as a criterion to 
obtain the 95% CL upper limit on 

In the bottom panel of Fig. we show, with a solid 
curve, the corresponding to the 95% CL upper limit 
for a = 2.2, 5 = 2, and z^ = 1.5, compared with the 
cross-correlation data between the > 1 GeV photons and 
the 2MASS galaxies, which gives the major contribution 
to the The dotted curve, in contrast, is the 95% 
CL upper limit due to the spectral data alone, and it is 
clearly inconsistent with the cross-correlation measure¬ 
ment. The top panel of the same figure shows the corre¬ 
sponding energy spectra for both approaches. It is clear 
that the source with the parameters adopted in Fig. [T] 
cannot be the main component of the IGRB spectrum be¬ 
cause the cross correlation provides a tighter constraint: 

£;95%CL ^ 5 9 ^ j^q 44 gj.g y^-1 Mpc-^. 

Figure]^ shows the 95% CL upper limits on as a 
function of a, for 5 = 2 and Zc = 1.5. For a wide range of 
spectral indices, the cross-correlation data provide con¬ 
straints more stringent by up to one order of magnitude 
than the spectral data. We also find that the difference 
is larger for smaller 5, since the cross correlation con- 



a 

FIG. 2. The 95% CL upper limits on the local gamma-ray 
luminosity density between 100 MeV and 100 GeV, as 
a function of a for 5 = 2 and Zc = 1.5. The limits due to 
spectrum and cross correlation data are shown as dotted and 
solid curves, respectively. 


straints are stronger for smaller redshifts, particularly 
due to the 2MASS galaxies. For <5 > 4, we find that both 
the spectrum and cross correlations provide comparable 
constraints on E-y^. The dependence on Zc, on the other 
hand, is significantly weaker as long as Zc > 1. 

We note that, to be conservative, we did not include 
secondary gamma rays that are generated by electromag¬ 
netic cascades, which would improve the spectral con¬ 
straints m- If the intergalactic magnetic fields are suffi¬ 
ciently weak such that the cascades do not produce halos 
or larger diffuse emission (e.g., m). the tomographic 
constraints will be also improved by the same factor. 

Constraints on high-energy neutrinos .—If neutrinos 
are produced by cosmic ray protons via pp interactions, 
their intensity is related to that of gamma rays [55] : 

« 6/^,no-EBL(L^7), (9) 

with Ey = 2 Ey. Here, Ii, is the neutrino intensity for all 
flavors, and Ty.no-EBL is the gamma-ray intensity without 
EBL absorption. Therefore, constraints on Ij{Ey) (or 
Ey^o), for each set of the parameters (a, 5, Zc), can be 
directly transformed into those of a neutrino intensity in 
the TeV-PeV energy range through Eq. (§. 

The top panel of Fig.[I]shows the 95% GL upper limits 
on I,j{Ey) from the IGRB spectrum (red dotted) and the 
cross-correlation data (red solid) for a = 2.2, 5 = 2, and 
Zc = 1.5. We find that while the IGRB spectrum analysis 
suggests this particular model to be compatible with the 
IceGube data, the tomographic approach constrains it as 
a subdominant source. 

Figure [^ shows the dependence of the neutrino inten¬ 
sity integrated above 25 TeV on a and 5 for a fixed 
Zc = 1.5. The intensity range preferred by the best-fit 
single power-law model of the IceGube data [5] is also 
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FIG. 3. The 95% CL upper limits on the neutrino intensity 
integrated above 25 TeV as a function of 5 for various values 
of a and fixed Zc = 1.5. Thick and thin curves show the limits 
due to the tomographic and spectral analyses of the IGRB, 
respectively. The horizontal magenta band shows the 68% CL 
interval of the best-ht single power-law model for the IceCube 
neutrino data [5], corresponding to the neutrino band shown 
in Fig. [g 


shown for comparison. For each model characterized by 
{a, 5), we show constraints due to the spectral and tomo¬ 
graphic data, as thin and thick curves, respectively. Note 
that the tomographic analysis gives tighter constraints by 
up to one order of magnitude with respect to the spec¬ 
tral analysis, especially for small S. In particular, for any 
source class slowly evolving (e.g., (5^3), even a very hard 
spectrum such as E~^'^ is nearly excluded as dominant 
source for the IceCube neutrinos. Any soft source with 
a > 2.2 should contribute much less to the total neutrino 
flux than previously expected (e.g.. Refs. [211147]). Mod¬ 
els with spectrum as hard as E~^, on the other hand, are 
still compatible with the IceCube flux level. 

Discussion and outlook .—Under the hypothesis that 
the TeV-PeV IceCube neutrinos are mostly generated 
from pp interactions in a single astrophysical source class 
(or more classes with similar properties). Fig. implies 
that a model with a « 2.15 and 5 « 4 (for Zc = 1.5) 
can explain most of the neutrino flux. At the same time, 
sources of this kind can explain most of the IGRB flux as 
well as the measured cross correlations. We note that in 
order for such a hard spectrum to be compatible with the 
IceCube data, a PeV spectral cutoff is required |S] (but 
data in the northern hemisphere still allow it without 
a cutoff i)- Otherwise, the comparison of the current 
data set with our results might suggest a mixed pp-pj, 
or even a pure py origin of the IceCube neutrino events. 

Interestingly, starburst galaxies well satisfy the above 


conditions for the pp origin, although efficient cosmic ray 
confinement needs to be achieved [191HH l2H] • While di¬ 
rect gamma-ray measurements of the redshift evolution of 
star-forming galaxies are not yet available, observations 
of their infrared luminosity (or of the star-formation rate) 
support such steep evolution. In particular, the evolu¬ 
tion of starburst galaxies is characterized by d > 4 up to 
Zc ^ 1.5 [58] . Here, we assumed that the local correlation 
between infrared and gamma-ray luminosities [63j holds 
also at high redshifts. 

Based on a modeling of resolved gamma-ray sources. 
Ref. [M] argued that about 20-30% of the IGRB above 
100 MeV can be explained by blazars (a subclass of 
AGNs). Furthermore, for energies above ~100 GeV, the 
blazar contribution can be substantial, explaining most 
of the IGRB data and leaving little room for any other 
source. This might point toward an even harder source 
population with steep redshift evolution for the neutri¬ 
nos, which would be, however, subdominant both in the 
IGRB flux and cross correlations. For example, in the 
case of a = 2 and ^ = 4, once we tune the gamma-ray lu¬ 
minosity density to match the level of ~10% of the IGRB 
flux and cross correlations, the same model could explain 
most of the neutrino data. 

Glusters and groups of galaxies have also been investi¬ 
gated as potential neutrino sources [301137] j where cosmic 
rays, generated through large-scale-structure shocks [371 
SO] or injected by star-forming galaxies [27], interact with 
the intracluster medium. Since the cluster/group num¬ 
ber density decreases as a function of redshift, imply¬ 
ing a small value of 6, tomographic constraints are very 
stringent. When considering starbursts or AGNs in clus¬ 
ters/groups, their quick redshift evolution has to be cou¬ 
pled with the negative one of clusters. As an example, 
we calculated that the overall evolution is locally charac¬ 
terized by (5 < 2 that quickly decreases to negative val¬ 
ues for z > 0.5. In addition, clusters are largely biased 
with respect to dark matter (i.e., 6-y ^ 5 for 10^^Mq 
and z — 0 [5S]), making the tomographic constraints 
tighter than those shown in Fig. ]^ Therefore, clusters 
and groups are disfavored by the cross-correlation data. 

These arguments cannot be applied to py sources, such 
as AGNs [T21[T3|[Tni[TnilIH]- This is because the threshold 
for py interactions is typically very high. It is also argued 
that such sources may be optically thick for GeV gamma 
rays [66]. In any case, it appears difficult that AGNs can 
be responsible for all the IceGube neutrino events. In 
fact. Ref. [T3] recently suggested that the diffuse emission 
from blazars can explain the IceGube neutrino flux at 
energies above ^PeV only. 

In conclusion, the tomographic method that we apply 
for the first time to high-energy neutrinos yields tight 
constraints on the properties of any hadronuclear source, 
providing complementary bounds on their injection spec¬ 
tral index and redshift evolution. In particular, we show 
that only hard spectrum sources with fast redshift evolu- 
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tion can produce a neutrino flux at the same level as the 
IceCube measurement. The potential relevance of this 
method in connection with high-energy neutrinos is ex¬ 
pected to quickly increase in the near future, because of 
the growing galaxy samples for the cross-correlation anal¬ 
ysis, including cosmic shear measurements that already 
seem promising [HThfiHj . 
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